# GMM in finance: a tutorial
# authors: Alan De Genaro
# created: 07 / Sept/2020


#####################################################################
# Set working directory 
this.dir <- dirname(parent.frame(2)$ofile)
setwd(this.dir)
######################################################################
# load dataset

Data_CCAPM.df = read.csv(file="Data_CCAPM.csv")
rownames(Data_CCAPM.df) = Data_CCAPM.df[,1]
Data_CCAPM.df = Data_CCAPM.df[,-1]

# Load gmm library
library(gmm)
options(digits=4, width=70)


############################################################################################
# Non-linear GMM estimation of equation asset pricing model
############################################################################################

# create excess return series for GMM estimation 
Data_CCAPM.mat = as.matrix(Data_CCAPM.df)
n.col = ncol(Data_CCAPM.mat)
excessRet.mat = apply(Data_CCAPM.mat[,2:(n.col-1)],2,
                      function(x,y){x-y},
                      Data_CCAPM.mat[,"RF"])

# data set for GMM estimation using J risky assets and a risk free asset
# data = (C(t+1)/C(t),1+Rf(t+1),R1(t+1)-Rf(t+1),...,
#         RJ(t+1)-Rf(t+1))
eulerDataJ.mat = cbind(Data_CCAPM.mat[,"CONS"],1+Data_CCAPM.mat[,"RF"],
                       excessRet.mat)
colnames(eulerDataJ.mat)[1] = "CONS"
colnames(eulerDataJ.mat)[2] = "RF"

# moment function for basic Euler equation asset pricing


eulerJ.moments <- function(parm,x=NULL) {
  # parm = (beta,gamma)
  # data = (C(t+1)/C(t),1+Rf(t+1),R1(t+1)-Rf(t+1),...,
  #         RJ(t+1)-Rf(t+1))
  n.col = ncol(x)
  sdf = parm[1]*x[,1]^(-parm[2])
  d1 = sdf*x[,2] - 1
  d2 = as.matrix(rep(sdf,(n.col-2))*x[,3:n.col])
  return(cbind(d1,d2))
}

start.vals = c(1,5)
names(start.vals) = c("beta","alpha")
euler.mom = eulerJ.moments(start.vals,eulerDataJ.mat)


#####################################################################
# gmm estimation
#####################################################################

# Two-Step GMM
start.vals = c(1,5)
names(start.vals) = c("beta","alpha")
eulerJ.gmm.fit2S = gmm(g=eulerJ.moments, x=eulerDataJ.mat, 
                     t0=start.vals, type="twoStep",
                     vcov="iid", optfct="optim")

# generates first column of "Tabela 1 (Table 1)" on paper
summary(eulerJ.gmm.fit2S)


# Iterated GMM
start.vals = c(1,5)
names(start.vals) = c("beta","alpha")
eulerJ.gmm.fit = gmm(g=eulerJ.moments, x=eulerDataJ.mat, 
                     t0=start.vals, type="iterative",
                     vcov="iid", optfct="optim")

# generates second column of "Tabela 1" on paper
summary(eulerJ.gmm.fit)

# CUE using iterated GMM as initial value
eulerJ.cugmm.fit = gmm(g=eulerJ.moments, x=eulerDataJ.mat, 
                     t0=coef(eulerJ.gmm.fit), type="cue",
                     vcov="iid", optfct="optim")

# generates third column of "Tabela 1" on paper
summary(eulerJ.cugmm.fit)

##########################################################################
# Including more instruments
#########################################################################

# moment function for basic Euler equation asset pricing

euler2J.moments <- function(parm,x=NULL) {
  # parm = (beta,gamma)
  # data = (C(t+1)/C(t),1+Rf(t+1),R1(t+1)-Rf(t+1),...,
  #         RJ(t+1)-Rf(t+1),C(t)/C(t-1))
  n.col = ncol(x)
  sdf = parm[1]*x[,1]^(-parm[2])
  d1 = sdf*x[,2] - 1
  d2 = as.matrix(rep(sdf,(n.col-3))*x[,3:(n.col-1)])
  d3 = d1*x[,n.col]
  d4 = d2*x[,n.col]
  return(cbind(d1,d2,d3,d4))
}

n.row = nrow(eulerDataJ.mat)
eulerData2J.mat = cbind(eulerDataJ.mat[2:n.row,], eulerDataJ.mat[1:(n.row-1), "CONS"])
colnames(eulerData2J.mat)[ncol(eulerData2J.mat)] = "CONS(-1)"

# check moments
start.vals = c(1,5)
names(start.vals) = c("beta","alpha")
euler.mom = euler2J.moments(start.vals,eulerData2J.mat)

#####################################################################
# gmm estimation
#####################################################################

# Two-Step GMM
start.vals = c(1,5)
names(start.vals) = c("beta","alpha")
euler2J.gmm.fit2S = gmm(g=euler2J.moments, x=eulerDataJ.mat, 
                       t0=start.vals, type="twoStep",
                       vcov="iid", optfct="optim")

# generates first column of "Tabela 2 (Table 2)" on paper
summary(euler2J.gmm.fit2S)

# Iterated GMM
start.vals = c(1,5)
names(start.vals) = c("beta","alpha")
euler2J.gmm.fit = gmm(g=euler2J.moments, x=eulerData2J.mat, 
                     t0=start.vals, type="iterative",
                     vcov="iid", optfct="optim")

# generates second column of "Tabela 2" on paper
summary(euler2J.gmm.fit)

# CUE using iterated GMM as initial value
euler2J.cugmm.fit = gmm(g=euler2J.moments, x=eulerData2J.mat, 
                       t0=coef(euler2J.gmm.fit), type="cue",
                       vcov="iid", optfct="optim")

# generates third column of "Tabela 2" on paper
summary(euler2J.cugmm.fit)

##################################################################

# Compute expected and predicted excess return.
# Use one-step GMM for the model with z={1} as instrument.

# compute average excess returns
excessRet.hat = colMeans(excessRet.mat)


# compute predicted excess returns from asset pricing model
predRet <- function(parm,x) {
   ncol = ncol(x)
  sdf = parm[1]*x[,1]^(-parm[2])
  tmp <- function(x,y) { -var(x,y)/mean(y) }
  ans <- apply(x[,3:ncol],2,FUN=tmp,sdf)
  return(ans) 
}


predRet.hat = predRet(eulerJ.gmm.fit2S$initTheta,eulerDataJ.mat)

# plot average excess returns against predicted returns
# annualized monthly data

predRet.hat.a = (1+predRet.hat)^12 - 1
excessRet.hat.a = (1+excessRet.hat)^12 - 1

# generates "Figura 1 (Figure 1)" on the paper 
plot(predRet.hat.a,excessRet.hat.a,ylab="Mean excess return",
     xlab="Predicted mean excess return",xlim=c(0,0.15),ylim=c(0,0.15), pch = 24,col="darkgreen", bg="darkgreen")
abline(a=0,b=1)


